Crack growth by surface diffusion in viscoelastic media 
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We discuss steady state crack growth in the spirit of a free boundary problem. It turns out that 
mode I and mode III situations are very different from each other: In particular, mode III exhibits a 
pronounced transition towards unstable crack growth at higher driving forces, and the behavior close 
to the Griffith point is determined entirely through crack surface dissipation, whereas in mode I the 
fracture energy is renormalized due to a remaining finite viscous dissipation. Intermediate mixed- 
mode scenarios allow steady state crack growth with higher velocities, leading to the conjecture that 
mode I cracks can be unstable with respect to a rotation of the crack front line. 
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The growth of cracks is major puzzle in solid state 
physics and materials science, still lacking a convincing 
physical description. It links macroscopic material prop- 
erties to microscopic effects in the tiny tip region and 
raises the important question which features of crack 
growth are generic and can be attributed to larger classes 
of materials. The latter has recently attracted inter- 
est through the phase field method, which reformulates 
crack propagation as a moving boundary problem, where 
the tip scale can be either determined through intrin- 
sic scales of the modeling technique or by macroscopic 
selection principles. Here, the Asaro-Tiller-Grinfeld in- 
stability (ATG) [J, i], which is usually understood as 
the morphological instability of a uniaxially stressed sur- 
face or interface, turns our to be strongly related to this 
problem: A perturbation of an initially straight surface 
increases the surface area, but diminishes the stored elas- 
tic energy, and therefore decreases the total energy of the 
system, provided that the wavelength of perturbation ex- 
ceeds a critical (macroscopic) lengthscale, which depends 
on the applied stress. We note that the same ingredients, 
a counterplay between a release of elastic energy and an 
increase of surface or fracture energy are the basic mech- 
anisms to understand the propagation of cracks beyond 
the Griffith point. Nevertheless, it is known that this in- 
stability leads to a breakdown of the physical description, 
as the unstable solid forms deep grooves, which, after a 
finite time, advance with infinitely high velocity and van- 
ishing tip radius. The reason for this breakdown is the 
absence of an additional microscopic lengthscale for selec- 
tion of a crack-like tip radius. Hence, understanding the 
selection of a crack tip radius in fracture has important 
implications also for the stability of stressed surfaces. 

One of the central questions for any crack model is 
the role of dissipation, which is directly connected to the 
quest for selection mechanisms for a tip scale. The elastic 
loading, which is applied far away from the crack tip, is 
usually only partially used to create the (macroscopically 
visible) crack surfaces; especially for higher propagation 



speeds a microbranching instability can significantly in- 
crease the fracture energy Q ■ This already indicates that 
the three-dimensional geometry is important for a full 
understanding of crack propagation, and it is one the 
goals of this paper to shed light on this by investigat- 
ing different modes of loading. In the tip region, the 
local temperature can rise significantly and even exceed 
the equilibrium melting temperature, which shows that 
dissipation can be very strong there. Although we have 
demonstrated that even pure dynamical linear elasticity 
can regularize the singular crack tip [J], it is natural to 
assume that deviations from a pure elastic behavior can 
play a crucial role, which can also contribute to dissi- 
pation; plasticity is an important example H, but still 
a full description has not yet been archived there. It is 
obvious that a very detailed modeling of the tip region 
is required to investigate these different effects, which in- 
clude a self-consistent selection of the crack shape itself. 

To address these important questions, we propose a de- 
scription of crack propagation in the spirit of interfacial 
pattern formation processes by inclusion of viscoelastic 
effects. This picture goes beyond the usual small scale 
yielding that is frequently used in the modeling of brittle 
fracture and includes two dissipative mechanisms: First, 
there is dissipation directly at the crack surface; the in- 
coming flow of elastic energy is partially converted to 
surface energy in order to advance the crack, and the re- 
maining part is converted to heat. Second, an extended 
zone of viscous dissipation is formed around the crack. 
We note that this problem is quite complicated as the 
shape of the crack, its velocity and the distribution be- 
tween viscous and interfacial dissipation have to be de- 
termined self-consistcntly. 

Viscous dissipation in mode I fracture has been dis- 
cussed in the literature, and although our results quali- 
tatively agree, our model makes a further step as it in- 
troduces this effect as way to intrinsically regularize the 
tip-singularity by selection of the crack tip radius. In con- 
trast, other models assume a Barenblatt crack tip model 
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or similar ad-hoc regularization criteria. For details see, 
for example, [|| and references therein. 

For simplicity we assume that the system obeys a 
translational invariance in one direction, thus it is ef- 
fectively two-dimensional. We assume an isotropic lin- 
ear viscoelastic medium, Ui and are displacement and 
strain respectively. The total stress, a ik = a^f' + <r^ ls ' , 
is decomposed into the elastic stress, which is given by 
Hooke's law (with elastic modulus E, Poisson ratio v), 
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which is related to the displacement rate through the 
viscosities r\ and Since we concentrate here on slow 
fracture with velocities far below the Raylcigh speed, the 
assumption of static viscoelasticity is legitimiate, thus 
daik/dxk — 0. On the crack contour, the total normal 



and shear stresses have to vanish, a r , 



0, with 



the interface normal and tangential directions n and s. 
The driving force for crack propagation is given by the 
chemical potential § 
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with 7 being the interfacial energy per unit area and f2 
the atomic volume; the interface curvature k is positive 
for a convex crack shape. Surface diffusion leads to the 
following expression for the normal velocity at each in- 
terface point 
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with the surface diffusion constant D (dimension m 4 /s). 
Notice that To := 2?y(l + v)/E defines a timescale, 
thus (Dto) 1 / 4 defines a lengthscale parameter which ul- 
timately leads to selection of the tip scale. In the gen- 
eral case, another timescale (which does not differ signif- 
icantly from To) is set similarly by the viscous coefficient 
£, and we discuss the specific case that these timescales 
are equal, i.e. C, = 2r;[u/(l — 2u) + 1/3]. Of course, this 
simplification is only relevant for mode I fracture, as the 
second scale does not appear in mode III. Altogether, the 
above set of equations fully defines the problem. 

We note that for steady state growth with velocity v, 
the last equation can be integrated once, and we obtain 
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We illustrate the procedure to solve the moving- 
boundary viscoelastic problem for mode III fracture in 



the steady state regime; for mode I loadings a similar 
approach can be used, which will be explained in de- 
tail elsewhere. The crack is located in the xy plane and 
propagates in positive x direction with velocity v. Then 
Newton's equation reads 



V 2 (u z + r u z ) = 0, 



(6) 



from which we obtain by differentiation ~S/ 2 <r xz = 
V 2 (Tj, z = 0. We represent the total stress through an 
analytical complex potential E with a xz — 3(E) and 

o yz = 3i(E). For steady state growth, a^ 1 ^ — — VTaa^ 3, 
and we therefore make a similar ansatz for the represen- 
tation of the elastic fields through an analytical function 
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This also 



guarantees the integrability of the strain field. The force 
balance Eq. ((6|) is then satisfied for solutions of the com- 
plex differential equation 
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with z — x + iy. For the total stress we use a multipole 
expansion with a branch cut along the negative real axis, 
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m — 1 



with real coefficients Ai. The main mode m = 1 is re- 
lated to the stress intensity factor, A\ — Km / ^,{2tt) 1 / 2 
(/i = £7/2(1 + v) is the shear modulus). The other coeffi- 
cients are adjusted such that the boundary condition on 
the (extended) crack shape, a nz — 0, is satisfied. To this 
end we minimize the residual stress functional J <J^ z ds 
(integrated along the crack contour) with respect to the 
expansion coefficients and solve the arising linear prob- 
lem numerically for a known crack shape. We restrict 
the calculation to a finite number of modes M in such 
a way that the final result does not change noticeably if 
the accuracy is increased. Eq. ([7]) can now be solved for 
each mode, and we obtain 
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with the recursion relation 
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, (10) 

(11) 



The integration constant is chosen such that far away 
from the crack tip the purely elastic behavior is retained. 
These expressions can then be used to obtain the strain 
from Hooke's law, and - after integration - the displace- 
ment field u,. 
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FIG. 1: Shape of a mode III crack for A = 2.5 in the steady 
state regime. The total incoming elastic energy flux is con- 
verted into surface energy, surface dissipation and viscous 
bulk dissipation, which is localized to the scale vto around 
the crack tip (visualized by the color coding). 



The strategy of solution is therefore as follows: First, 
for known crack shape, the total stress problem is solved 
in the spirit of Eq. (jSJ), delivering the coefficients of ex- 
pansion A4. They, together with a given crack speed 
are used to determine the elastic stress field according to 
Eqs. ([§))- (fTTj) . Then, in the next step, the chemical po- 
tential © can be computed using Hooke's law. Finally, 
the steady state equation © is a nonlocal and nonlinear 
relation which is used to determine a new guess for the 
crack shape and velocity. With them, the whole proce- 
dure is iterated until a self-consistent solution is found. 

We define a dimensionless driving force 



A = A, 
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where we already included the possibility of mixed-mode 
loading, and A = 1 is the Griffith point. From now on, 
we set v = 1/3. Fig. [1] shows a typical steady state crack 
shape for mode III loading in the reference frame (La- 
grangian coordinates), i.e. the elastic displacement is not 
included. First, we clearly see that the crack tip scale is 
selected self-consistently, and the finite time cusp singu- 
larity of the ATG instability does not occur. Therefore, 
the presence of viscous bulk dissipation is a way to cure 
this well-known problem. Second, it is important that far 
behind the crack tip the opening decays to zero, which is 
a consequence of mass conservation, as expressed by the 
equation of motion for surface diffusion Diffusive 
transport is therefore restricted to the tip region, and no 
long-range transport is required. Qualitatively, the crack 
shapes for mode I look very similar. 

Interestingly, the propagation velocity differs quite sig- 
nificantly for mode I and mode III fracture, as shown in 
Fig. [2j For mode III, the crack speed increases with the 
driving force, until it reaches a maximum at A w 3.5, 
then it decreases, and obviously steady state solutions 
do not exist beyond the point A w 3.8, where the stable 
branch merges with another (unstable) solution. Beyond 



0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




FIG. 2: Steady state propagation velocity as function of the 
driving force for pure mode III and mode I fracture. Addition- 
ally, mixed mode situations with A//A = 0.15 and Aj/A = 
0.85 are displayed. The velocity scale is «o = (Dr ' 3 ) 1 ^ 4 ' . 
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FIG. 3: Half crack height as function of the driving force. The 
lengthscale used here is ho = (Dto) 1//4 . 



the bifurcation point we expect crack branching, in anal- 
ogy to our findings for fast brittle fracture [4[. 

Fig. [3] shows the maximum height of the crack as func- 
tion of the driving force for different loadings. At A « 1.1 
the size of the mode III steady state crack diverges and 
v — ► 0. The viscous dissipation becomes negligible here, 
but the surface dissipation remains finite. This point can 
be interpreted as the point of ductile-to-brittle transition: 
Below it the size grows indefinitely in time and the crack 
slows down, while above this point steady state solutions 
with a finite tip scale exist. 

Starting from a pure mode III crack, we can now in- 
clude additional mode I loadings. Fig. [2] shows that this 
shifts the bifurcation point towards higher values and 
therefore extends the range of steady state solutions to- 
wards higher driving forces. Again, the crack blunts close 
to the 'nominal' Griffith point A = 1. Simultaneously, 
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FIG. 4: Distribution of energy consumption. The fraction 
above the solid curve is the relative contribution of surface 
energy generation, the part below the dashed line the viscous 
dissipation (shown here for two different admixtures of mode 
I loading; for 15% mode I also for the unstable branch). The 
remaining part between the curves is the surface dissipation. 



the propagation velocity is significantly reduced in the 
regime of small A, as can be clearly seen in the compari- 
son between the cases with 15% and 85% mode I contri- 
bution. Effectively, this establishes an interval of driving 
forces, where the crack speed is very low, and only af- 
ter this plateau it sharply increases; this effect becomes 
more pronounced as the crack loading is more mode I 
dominated. The same plateau can also be found in the 
tip scale, see Fig. [3] 

For the case of mode I, finally, steady state solutions do 
not exist below A = 2.6; this result has to be interpreted 
as a limiting case with very slow creep with velocities and 
tip radii significantly lower than above the point A = 2.6. 
Literally, of course, growth starts at A = 1 due to energy 
conservation. The presence of this plateau is quite re- 
markable, as this effectively renormalizes the 'apparent' 
Griffith point - the driving force where the velocity starts 
to increase sharply - to a substantially higher value than 
A = 1, although the viscous dissipation remains finite 
on the 'creep branch'. Again, the crack speed increases 
monotonically with the driving force, and the bifurca- 
tion to unstable growth occurs only at very high driving 
forces. 

Finally, Fig. [4] shows how the different mechanisms 
contribute to the total energy consumption: From the 
total dimensionless driving force the amount A — 1 is 
dissipated, and the remaining part is used for the cre- 
ation of the surfaces of the advancing crack. Obviously, 
this contribution becomes less important in comparison 
to the true dissipation for higher driving forces. The vis- 
cous dissipation is 

l^ = ±f RdV=±f a<Z>iH*d. (13) 



with the Rayleigh dissipation function R = e^, and 
the integration domain V is the solid phase. The latter 
equality in Eq. (fT3|) holds in steady state, where S is the 
crack contour. The surface dissipation is 

A s = J o [ ^\ lk n x ds - 1 = ^ J y 2 ds, (14) 

with the horizontal component of the interface normal 
n x ; again, the latter expression, which follows from 
Eq. ([5]), is valid only in the steady state regime. Al- 
together, we have A = A s + A.„ + 1. 

Starting from the Griffith point the viscous dissipation 
continuously increases up to the point where the stable 
steady state solution branch terminates. It is quite re- 
markable, that for mode I dominated cracks the viscous 
dissipation A v is much larger than A s , which shows that 
bulk dissipation can indeed play a crucial role. Notice 
that these (dimensionless) predictions do not depend on 
model parameters. 

The obtained results lead to the striking conclusion, 
that the apparent Griffith point may depend quite sub- 
stantially on the mode of loading. Although most mod- 
els in the literature are discussed either in the mode I or 
mode III case only, we clearly see here that the behavior 
can be significantly different in these cases, as soon as 
bulk dissipation is taken into account. For the specific 
case of crack propagation in viscoelastic media we obtain 
that the onset of steady state growth is shifted towards 
higher values in mode I. This leads to the interesting con- 
sequence that by a rotation of the crack front, which can 
induce a mode III stress intensity factor, the apparent 
Griffith threshold can be reduced and the crack speed 
increased; a fully time-dependent three-dimensional sim- 
ulation should shed light on this conjecture. 

In summary, we developed a model for crack propa- 
gation in viscoelastic media in the spirit of an interfa- 
cial pattern formation process. Motion occurs due to 
surface diffusion along the extended crack shape, which 
is - together with the propagation velocity and the tip 
scale - selected self-consistently. The steady state regime 
of the model is solved numerically using a series expan- 
sion method and a sharp interface description, which effi- 
ciently separates the microscopic crack tip scale from the 
system size. The results show that the bulk dissipation 
in the surrounding of the crack tip can play a substantial 
role especially for higher driving forces, and the crack 
velocity depends crucially on the mode of loading. 
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